%this version of main aims to incorporate features retrieval

%created 9/8/20 as mainv02 in what is now v01 to incorporate
%features retrieval in Lehman example
%10/??/20 features retrieval incorporated into Depression Babies
%10/21/20 made Lehman notation consistent with new notation in
%text, added figures to respond to Referee 5, edited figures for clarity

clear all

close all

zeta = 0.35

%options
FS = 18;
featuresretrieval = 1


model = 2;







if model == 1
    
    tau_exp = 4; %experience interval for model 1.
    
    load paverage
    
    p = .02 %probability of a financial disaster
            %tau = 0;
    tau = 1;
    tvector = [0:1:49];  %posterior sample size 
    
    p_prior = .5 %the prior value of p
    pbar = 1/2; % has the interpretation of an upper limit
                %on p, used to benchmark y and sigma
    

    alpha = p_prior*tau+1  %observed disasters in prior sample
    beta = (1-p_prior)*tau + 1;
    
    
    
    %posterior mean of p
    p_post = (p*tvector + alpha)./(tvector + alpha +beta);
    
    %for zeta<1
    pzeta = zeros(1,length(tvector));
    pneglect= pzeta;
    
    pzeta(1) = p_prior;
    pneglect(1) = p_prior;
    for t = 2:length(tvector)
        pneglect(t)  = pneglect(t-1)*(1-zeta);
    end

    if featuresretrieval
        %this assumes that we have already run
        %maindepressionbabies!  Make sure parameters are the same.
        pzeta = paverage;
    else
        for t = 2:length(tvector)
            pzeta(t) = pzeta(t-1)*(1-(1-zeta)/(tau_exp*t));
        end
    end
    
    
    
    %portfolio allocation
    y = 1/(1-pbar);
    sigma = sqrt(pbar/(1-pbar));

    %Bayesian allocation
    thetapost = (1-p_post*y*sigma)/sigma^2;
    %FI allocation
    thetapostcorrect = (1-p*y*sigma)/sigma^2;
    %RCT -- prior does not change.
    thetarcont = (1-p_prior*y*sigma)/sigma^2;
    
    thetazeta = (1-pzeta*y*sigma)/sigma^2;
    thetaneglect = (1-pneglect*y*sigma)/sigma^2;
    
    
     figure(1)
     %h = plot(tvector,p_post,'--',tvector,p_prior*ones(size(tvector)),'-',tvector,pzeta,'-o',tvector,pneglect,'-x',tvector,paverage,'+');
     h = plot(tvector,p_prior*ones(size(tvector)),'-o',tvector,pzeta,'-x',tvector,p_post,'--',tvector,pneglect,'-');
    set(gca,'fontsize',FS)
    set(h,'linewidth',2)
    axis([-1 50 0 .52])
   
    ylabel('posterior probability','fontsize',FS)
    xlabel('years of data','fontsize', FS)
    legend('Retrieved context (\zeta = 1)',...
           'Retrieved context (\zeta = 0.35)','Bayes','Neglected risk');
    %print -depsc posteriorprobability.eps
    print -deps posteriorprobability.eps  

    
    figure(2)
    h = plot(tvector,  thetarcont*ones(size(tvector)),'-o', ...
         tvector, thetazeta,'-x',tvector,thetapost, '--',tvector,thetaneglect,'-');
    set(gca,'fontsize',FS)
    set(h,'linewidth',2)
        axis([-1 50 -.02 1.02])
    ylabel('allocation','fontsize',FS)
    xlabel('years of data','fontsize', FS)
    %print -depsc portfolioallocation.eps
    print -deps portfolioallocation.eps
    
    
end


if model == 2
   
    
    %lehman
    pc = .05  %probability of a crisis
    q = .5

   
    tau = 100; %length of prior sample
    t0 = 3;  %periods of crisis (September, October, November)
    t1 = 5;  %periods of recovery
    tprev = 2; %periods prior to crisis
 
    
    %utility parameters are standard
    beta = exp(-.02)
    etabar = -log(.85)
    muC = 0;
    muD = 0;
    lambda = 2
   
    %feature vectors
    f_series = [0;1;0]*ones(1,t0);  %crisis
    f_series = [f_series [1;0;0]*ones(1,t1)]; %post-crisis
    t = size(f_series,2);
    
    M0 = [1-pc,0,0; 0, pc*(1-q),pc*q];
    c0 = [1;0]
    m = size(c0,1)  %number of orthogonal contexts
    n = size(M0,2)  %number of features
    
    %initialize
    M_series = zeros(m,n,t+1);
    c_series = zeros(m,t+1);
    fin_series = zeros(n,t);
    PD_series = zeros(1,t);  %append PD0 later
    r_series = zeros(1,t);
    ErS_series = zeros(1,t);
    p_series = zeros(1,t);
    fins = zeros(n,m);
    
    %features retrieval step
    fin0 = [1;0;0];
    p0 = fin0(3);
    
    %inital values
    [PD0,r0,ErS0]  = PDdisasters(p0,beta,etabar,muC,muD,lambda)
    [PDfull,rfull,Erfull] = PDdisasters(pc*q,beta,etabar,muC,muD,lambda)
     
    c_series(:,1) = c0;  %first context will be time 0
    M_series(:,:,1) = M0;    
   
  
    id = eye(m);  %identity matrix the length of the
                           %context vector for use in creating
                           %retrieved features.
    cfsum = 0;
    for s= 1:t
        s = s
        Ms = M_series(:,:,s);
         cin = Ms*f_series(:,s);
        cin = cin/sum(cin);
        c_series(:,s+1) = (1-zeta)*c_series(:,s) + zeta*cin;
        cs = c_series(:,s+1);
        
        %features retrieval step
        for k = 1:m
            ek = id(:,k);  %kth basis vector
            fins(:,k) = (Ms'*ek)/sum(Ms'*ek);
        end
        finsum = fins*cs;
        fin_series(:,s) = finsum;
        p = finsum(3)
        [PD_series(s),r_series(s),ErS_series(s)] = ...
            PDdisasters(p,beta,etabar,muC,muD,lambda);

        %update memory, for now assume we use retrieved features
        %important, c_series and M_series include the initial
        %condition in their first entry.
        cfsum = cfsum + c_series(:,s+1)*fin_series(:,s)';
        M_series(:,:,s+1) = tau/(s+tau)*M0 + 1/(s+tau)*cfsum;
        p_series(s) = p;
    end
    
    PD_series = [PD0*ones(1,tprev),PD_series];
    r_series_unc = [r0*ones(1,tprev), r_series];
    r_series = max(r_series_unc,0);
    ErS_series = [ErS0*ones(1,tprev),ErS_series];
    p_series = [p0*ones(1,tprev),p_series];
   
    %    [PD_Bayes, r_Bayes] =  PDdisasters(pc*qtrue,beta,gamma,mu,sigma,xi,lambda)
    
    %note, no actual disasters
    R_series = exp(muD)*(PD_series(2:end)+1)./PD_series(1:end-1)-1;

    
    figure(1)
     
   
    yyaxis left 
    h = plot([-tprev+1:1:t],PD_series, [-tprev+1:1:t], ones(size([-tprev+1:1:t]))*PDfull,'--');
    xlabel('time','fontsize',FS)
    ylabel('Price-dividend ratio','fontsize',FS)
    set(gca,'fontsize',FS)
    axis([-inf, t1, min(PD_series)-3, max(PD_series)+1])
    set(h,'linewidth',2)
    %    print -deps PDLehman.eps
    
 
    yyaxis right
     h = plot([-tprev+1:1:t], r_series,'-',[-tprev+1:1:t],ones(size([-tprev+1:1:t]))*rfull,'--',[-tprev+1:1:t],zeros(size([-tprev+1:1:t])),':');
    xlabel('time','fontsize',FS)
    ylabel('riskfree rate','fontsize',FS)
    %legend('RCT','FI')
    set(gca,'fontsize',FS)
     axis([-inf, inf, min(r_series-.01), max(r_series+.04)])
        set(h,'linewidth',2)
    print -deps LehmanPD.eps
    

    figure(2)
    plot([-tprev+2:1:t],R_series,'o','markersize',10, ...
             'markerfacecolor','blue')
    hold on
    plot([-tprev+1:1:t],zeros(size([-tprev+1:1:t])),':','linewidth',2); 
    xlabel('time','fontsize',FS);
    ylabel('stock return','fontsize',FS);
    set(gca,'fontsize',FS)
    axis([-tprev+1, inf, -inf, inf])
    %[-tprev+2:1:t], ...
    %         r_series(1:end-1),'.-',[-tprev+2:1:t], ErS_series(1:end-1));
    %set(h,'markersize',10);
    %    set(h,'MarkerFaceColor',[1 .6. .6]);
    
    %test = p_series*(-exp(-(lambda-1)*etabar) + exp(-lambda*etabar) ...
    %                 + exp(etabar)-1)
    %test2 = ...
    %    -log(1+p_series*(exp(-(lambda-1)*etabar)-1))...
    %      + log(1+ p_series*(exp(-lambda*etabar)-1)) ...
    %            + log(1+ p_series*(exp(etabar)-1))
                                                    
    print -deps LehmanRS.eps
    
end


if model ==3 
    %this is the set-up for momentum
    p = .5; %probability of a high growth firm
    gH = .01
    gL = 0;
    %r = .10;
    r = .05
    T = 10;
    %initialize context vectors, prices and returns
    cW = zeros(1,T);
    cL = cW;
   
    PDW = cW;
    PDL = cW;
    
    %initial context
    c0 = [p,1-p]';
    %initial PD ratio
    PD0 = PD(c0(1),gH,gL,r)
    
    %cin for winners
    cin1W = [1,0]';
    %c1 winners
    c1W = rho*c0 + zeta*cin1W
    PD1W = PD(c1W(1),gH,gL,r)
    PricechangeW = (1+gH)*PD1W/PD0
       
    cin1L = [0,1]';
    c1L = rho*c0 + zeta*cin1L
    PD1L = PD(c1L(1),gH,gL,r)
    PricechangeL = (1+gL)*PD1L/PD0
    
    cW(1) = c1W(1);
    cL(1) = c1L(1);
    
    for t= 2:T
         cW(t) = rho*cW(t-1) + zeta;
         cL(t) = rho*cL(t-1);
    end
    

    cW = [c0(1),cW];
    cL = [c0(1),cL];
    PDW = PD(cW,gH,gL,r);
    PDL = PD(cL,gH,gL,r);
    RW = (PDW(2:end)+1)./PDW(1:end-1)*(1+gH)
    RL = (PDL(2:end)+1)./PDL(1:end-1)*(1+gL)
    
    PDWBayes = PD(1,gH,gL,r);
    PDLBayes = PD(0,gH,gL,r);
    
    tvector = [0:1:T];
    figure(1)
    h = plot(tvector,PDW,'-',tvector, PDL,'-', tvector, PDWBayes*ones(1,T+1), ':',...
         tvector,PDLBayes*ones(1,T+1),':');
    set(gca,'fontsize',FS)
    xlabel('time','fontsize',FS)
    axis([-1 T 19.5 26])
    gtext('Winners','fontsize',FS)
    gtext('Losers','fontsize',FS)
    set(h,'linewidth',2)
    print -deps momentumPD.eps
   
    figure(2)
    h = plot(tvector(2:end),RW-1,'-', tvector(2:end),r* ...
         ones(1,T),':', tvector(2:end),RL-1, '-');
    
        set(gca,'fontsize',FS)
    xlabel('time','fontsize',FS)
     legend('Retrieved context', 'Full information')
     gtext('Winners','fontsize',FS)
     gtext('Losers','fontsize',FS)
      set(h,'linewidth',2)
         print -deps momentumR.eps

    %now do the return for the second period, but generalize to
    %n periods.
    
    % cin2W = [1,0]';
    % c2W = rho*c1W + zeta*cin2W
    % PD2W = PD(c2W,gH,gL,r)
    % RW = ((PD2W+1)/PD1W)*(1+gH)
    
    % cin2L = [0,1]';
    % c2L = rho*c1L + zeta*cin2L
    % PD2L = PD(c2L,gH,gL,r)
    % RL = ((PD2L+1)/PD1L)*(1+gL)
    

    
end


if model == 4
    %this is the set-up for scary movies
    p = .02;  %probability of a risk
    q = .5;

    %values for y
    y_normal = 0;
    y_disaster = -.8;
    mu = .04;
    sig = .2;
    q_r= .5;
    
     r_1 = mu + sig;
    r_2 = mu - sig;
   
    P = [1-p 0 0; 0 p*(1-q) p*q];
    
    % c0 = P*[1 0 0]'

    % c0 = c0/sum(c0)

    c0 = [1-p p]';

    cin = P*[0,1,0]'
    cin = cin/sum(cin)

    c1 = (1-zeta)*c0 + zeta*cin
    
    thetas = [0:.05:1]
    V1 = zeros(size(thetas));
    V0 = V1;
    Vcorrect = V1;

    p = c0(2)
     for j = 1:length(thetas)
        theta = thetas(j);
        V0(j) = scarymovieu(theta,p,y_normal,y_disaster,q_r,r_1,r_2);
    end
    
    p = c1(2)
    for j = 1:length(thetas)
        theta = thetas(j);
        V1(j) = scarymovieu(theta,p,y_normal,y_disaster,q_r,r_1,r_2);
    end
    
 
    

    figure(1)
    h1 = plot(thetas,V0);
    a1 = gca;
    xlabel('Allocation','fontsize',FS)
    ylabel('Expected utility','fontsize',FS)
    figure(2)
    h2 = plot(thetas, V1);
    a2 = gca;
    xlabel('Allocation','fontsize',FS)
    ylabel('Expected utility','fontsize',FS)

    
    set([a1 a2],'fontsize',FS)
    figure(1)
    print -deps utility0.eps
    figure(2)
    print -deps utility1.eps
  
    

end



